We will use data from the article Determinants of healthcare worker turnover in intensive care units: A micro-macro multilevel analysis by Daouda, Hocine & Temime (2021).

Variables:

• stress level … DV (nurse’s reported stress level)
• support from colleagues … L1 IV (per nurse)
• type of ICU … L2 IV (per unit)

To disentangle L1 and L2 effect of support from colleagues we will group mean center it and also use its group means as a predictor.

So, instead of support from colleagues we will use:

• support_l1
• support_l2

Further, support_l2 will be grand mean centered.

• start simple (establish a baseline for evaluating larger models)
• add covariates one at a time, level by level, maybe centering certain variables
• finally, examine the random effects and variance components, beginning with a full set of error terms and then removing covariance terms and variance terms where advisable (for instance, when parameter estimates are failing to converge or producing impossible or unlikely values)

This strategy follows closely Raudenbush & Bryk (2002) approach.

Other strategies are possible. E.g., Diggle et al. (2002) begins with a saturated fixed effects model, determines variance components based on that, and then simplifies the fixed part of the model after fixing the random part.

fixed effects

• levels of a factor we want to draw inferences from
• would not change in replications of the study
• e.g., type of ICU

random effects

• levels of a factor which is just a sample from a larger population of factor levels
• not interested in drawing conclusions about specific levels
• but, interested in accounting for the influence of the random effect
• e.g., specific unit

Getting p-values in multilevel models is not trivial, because the exact distribution of the test statistics under the null hypothesis (no fixed effect) is unknown (Bates et al. 2015).

Approaches:

• t-values (ratios of parameter estimates to estimated standard errors) with absolute value above 2 indicate significant evidence that a particular model parameter is different than 0
• tests based on conservative assumptions, large-sample results, or approximate degrees of freedom for a t-distribution
• bootstrap, e.g.,
  • parametric (Efron, 2012)
  • residual (Carpenter, Goldstein and Rasbash, 2003)
  • wild (Modugno and Giannerini, 2015)

model_3 <- lmer(`stress level` ~ `type of ICU` + support_l1 + support_l2 +
  (support_l1 | CodeService), REML=T, data=df_ICU)

model_3_ML <- lmer(`stress level` ~ `type of ICU` + support_l1 + support_l2 +
  (support_l1 | CodeService), REML=F, data=df_ICU)

##  Groups      Name        Variance Std.Dev. Corr  
##  CodeService (Intercept)  0.59081 0.76864        
##              support_l1   0.19536 0.44200  -0.640
##  Residual                21.21546 4.60602

##  Number of Level Two groups =  30

##                   Estimate Std. Error     t value
## (Intercept)    22.72183488  0.4168353 54.51034544
## `type of ICU`1 -0.65742457  0.5796305 -1.13421329
## `type of ICU`2  0.93917067  0.6860464  1.36896077
## support_l1     -0.01408992  0.1409538 -0.09996129
## support_l2     -0.13764067  0.3798257 -0.36237851

Neither, ICU Type 1 or 2 seem to differ significantly from Type 0. But it seems like 1 and 2 could differ.

How would we test that?

Let’s try bootstrapped likelihood ratio test (BLRT).

b3 <- lmeresampler::bootstrap(model_3_ML, .f = logLik, type = "parametric", B = 1000)
b4 <- lmeresampler::bootstrap(model_4_ML, .f = logLik, type = "parametric", B = 1000)
lrt_b <- -2 * b4$replicates + 2 * b3$replicates
quantile(lrt_b, probs = c(.025, .975))

##      2.5%     97.5% 
## -81.85642  98.24528

CI does include 0 -> no model is significantly better.

-> ICs indicate better fit of model_4 and, although more parsimonious, it doesn’t fit the data significantly worse than model_3.

If your theory doesn’t emphasize the random slope, than just the random intercept seems to be sufficient.

Always pick the model that fits your theory.

Note: Theory is fixed before the analysis.